Mode coupling theory for sheared granular liquids 
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Abstract. Sheared granular liquids are studied by the mode coupling theory. It is shown that, in contrast to thermostatted 
systems, current correlations play an essential role in the dynamics. The theory predicts that the plateau of the density time- 
correlator disappears for most situations, while it appears in the elastic limit. The result is compatible with molecular dynamics 
simulations. 
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INTRODUCTION 

It is believed in the description of granular flows that 
the kinetic theory of Boltzmann-Enskog is applicable up 
to volume fraction <p < 0.5 [1, 2, 3, 4, 5]. However, any 
reliable thoery for dense granular flows with <p > 0.5 
does not exist, where the picture of instantaneous colli- 
sion breaks down. Experiments and simulations suggest 
that the constitutive equation of the shear stress exhibits 
a crossover from the Bagnold law to a power-law behav- 
ior with a non-trivial exponent, and eventually results in 
a yield stress for large <p [6, 7]. To understand this behav- 
ior by a first-principle theory is desired. 

On the other hand, Liu and Nagel [8] suggested the 
existence of a closed relationship between the jamming 
transition and the glass transition. Since then, to clarify 
this relation has been one of the hottest subjects in both 
granular and glassy physics [9]. This situation has moti- 
vated us to apply the mode coupling theory (MCT) [10] 
successfully used for glassy materials to granular materi- 
als [11, 12]. This is in accordance with the establishment 
of the notion of "granular liquids" [13, 14, 15, 16], where 
long-time, long-range correlations are significant, as well 
as in molecular liquids. Hayakawa and Otsuki [11] intro- 
duced a MCT for sheared granular liquids of hard-core 
inelastic spheres, but realistic granular grains are soft and 
the application of the pseudo-Liouvillian [1] is unneces- 
sary. Kranz et al. [12] studied MCT for driven inelastic 
spheres with a Gaussian white-noise thermostat, but the 
correspondence of their model with the actual vibrating 
granular systems, which incorporate non-Gaussian me- 
chanical forces, is quite unclear. 

With the above situation in consideration, we attempt 
to construct a liquid theory for sheared granular spheres 
with a soft-core potential, starting from the Liouville 
equation [17, 18, 19, 20]. In particular, we apply the 
MCT to obtain a set of closed equations for the time- 



correlators. Although MCT is still not completely estab- 
lished as a theory of glasses, its problem is not serious at 
least for sheared systems, so we adopt the conventional 
approach of projection operators and the Mori-type equa- 
tions [21, 22]. The physical significance of the current 
correlation is demonstrated by exhibiting its effect on the 
slow dynamics. 

FORMULATION 

In this section, we present the equation of motions for 
the time-correlators (MCT equation). Detailed deriva- 
tions for the basic formulation of sheared underdamped 
MCT can be found in Ref. [23], so we focus on the issues 
specific to granular systems. 

Our microscopic starting point is the SLLOD equation 
[24] for a system under a uniform and stationary shear 
with N spherical particles of mass m and diameter d, 

hit) = Z®- + K>r t (t), (1) 
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P,(t) = FfHt)+Ff s \t)-K- Pl (t), (2) 

where K"^ = y8^ x 8^ y is the XpL component of the shear- 
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rate tensor {X,jX — x,y,z), j is the shear rate, F \ is 

the conservative force with a soft-core potential, F| dls ' = 

Y,j^i^if^ is tne viscous dissipative force, where F^ 1 ^ 
is the pairwise contact disspative force, 
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with Vij ee Vi - v/, nj ee n - r h nj = |r l7 |, fij ee r u /r ih 
and &(x) is the Heaviside's step function, where &(x) = 
1 for x > and @(x) = otherwise. Here, £ is the bare 



viscous coefficient which is directly related to the restitu- 
tion coefficient. We consider frictionless granular materi- 
als to avoid the complexity to treat the Coulomb friction 
and the rotation of grains. The uniform shear velocity (in 
the x-direction) is v s h(y) = yy, and the distance of the 
two shear boundaries is denoted as L. Hence, the veloc- 
ity at the boundaries may satisfy = ±Lf/2. The sys- 
tem is assumed to be in equilibrium before time t = 0, 
and shearing and interparticle dissipation are turned on 
at t = 0, after which a steady-state is asymptotically 
reached for t — > °°. 

It is possible to recast the SLLOD equation, Eqs. (1) 
and (2), into the Liouville equation, and successively by 
the application of the projection operator [20], £?(t) = 

(*«*(,)) / ( NS m) n m + E* / (M) j x m 

to derive the Mori-type equation for the time-correlators. 
Here, n q = E£i«* r ' ~NS qfi and j\ = Zf =l pf e** /m 
are the density and current-density fluctuations, re spec - 
tively, with q(t) = (q x ,q y - jtq x ,q z ), and v T = \/T/m 
is the thermal velocity at the initial equilibrium state. In 
granular liquids, it is essential to include the projection 
onto the density-current modes £P n j(t), where 
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to the second projection operator ^(t) as [25] 

& 2 (t) = &nn{t) + &nj(t) 



(5) 
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to obtain a closure for the Mori-type equation, where 
^nn(t)X = Lk> P (Xn* m n p{t) ) I (N 2 S k{t) S p{t) ) n k(t) n p[t) 
is the conventional projection onto the pair-density 
modes. This choice of requires us to consider, 

in addition to the conventional time-correlators <t> q (t ) = 

(n q (t){t)n* q ) /N and H\{t 
tional time-correlators, 
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These time-correlators turn out to play essential and 
significant roles in granular liquids. 



Isotropic approximation 

It is almost hopeless to deal with four time-correlators 
with tensor indices. Hence, we attempt to reduce the 
degrees of freedom by applying the isotropic approx- 
imation. Since dissipation occurs irrespective of the 



anisotropy of the strain, application of this approxima- 
tion is expected not to be fatal. 

Based on the isotropic approximation, it is able to 
reduce the vector/tensor time-correlators into two scalar 
time-correlators, {4> (? (f), l I , (? (f)}. In particular, H^(t) is 

reduced to H%(t) ~ q l (t)j- t <S> q (t) / q(t) 2 as usual [23], 

andff^(f) to [26] 

fl^(/)~VW%W. (9) 

The correlator C q ^(t) is reduced to C q ^(t) ~ C g (t)8 1 ^ 
[26], where C q (t) is written in terms of as 
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The resulting MCT equations in the weak-shear regime 
read 
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■ f ' drr 2 g(r)&{r)[\ - j (qr) +2j 2 (qr)} , (13) 
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with the spherical Bessel function j n (x) (n = 0,2), and 
A XX = 4n n r drr 2 g{r) ^ {r) [j _ jQ ( qr)] (14) 

are the effective viscous coefficients. Here, n is the aver- 
age density. In the weak-shear regime, we keep terms up 
to linear order in the shear rate y and the bare viscous co- 
efficient £, and neglect quadratic and higher terms. For 
instance, C q (t) is approximated as C q (t) ~ <F¥ q {t)/dt. 
Note that J^(r) incorporates £, cf. Eq. (4). According to 
the choice of 3?2{t) in Eq. (6), there appear four terms in 
the memory kernels M q (x) and M q ^(x), respectively. A 
complete description of these terms is rather lengthy, so 
we report it elsewhere [25], and only show the explicit 
forms of M q {x) for illustration. Among the four terms in 
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FIGURE 1. Numerical results for the density time-correlator (left) and the density-current time-correlator (right). 



the first term M q (t) is the conventional term quadratic 
in* ? (f), 

which is responsible for the plateau of the density time- 
correlator. The other three terms originate in dissipation, 
where the two of them are given by 

« [T) ~ q 2 J (271)3 p ( T )2 V).*M.P(T)V*,P 

X* t (T)^*p(T), (17) 

M (3) (t) ~ -1 / -i!ly( eI ) „A y (vis)A 

x<&*(t)^(t), (18) 
and the remaining Mq{%) is neglected since it is 

(el) 

quadratic in £. Here, = (q ■ k)ct + (flf • p)c p , with 

p = q — k and the direct correlation function q, is the 
conventional vertex function, and V^™ p is the dissipative 
vertex function, whose explicit form is given by 

VS? = lffrg(rWr)*£S (>'-*<*') . (19) 

Note that Eq. (19) includes g(r),^(r) in the integrand, 
similarly to Eqs. (13) and (14). An important property 
of the dissipative memory kernels is that Mq (t) and 

M^q\x) are of opposite sign with respect to Mq {%). 
This implies that the interparticle dissipation operates to 
destroy the plateau of the density time-correlator in the 
j3 -relaxation regime. These features are commonly seen 
in the three terms of M^ (t) as well. 



NUMERICAL CALCULATION 

Now we present the result of the numerical calculation. 
We first describe the calculational conditions, starting 
with the non-dimensionalization scheme. The units of 
mass, length, and time are chosen to be m, d, and To = 

d/\>o, respectively, where vo = vjj — vjj ' = Ly is the 
relative velocity between the shear boundaries. In these 
units, the shear rate f= vq/L is non-dimensionalized as 
7* = 71b = d/L, which we require to be small enough, 
i.e. d/L <C 1 [27]. In the remainder, non-dimensionalized 
quantities are denoted with superscript * . 

Next, it is necessary to establish the relation between 
the bare viscous coefficient £ and the repulsive coeffi- 
cient e to compare our result with molecular dynamics 
simulations. As explained in Ref. [28], e is related to 
C* as e = exp [-£*T*], where x* = % / ' v '2k* - £* 2 is the 
contact duration time of spheres, if the contact force can 
be approximated by the linear spring model. Here, k* is 
the elastic spring coefficient of F' e1 '. In granular liquids, 
the steady-state temperature is determined by the balance 
of the work by shearing and the interparticle dissipation. 
The case of interest might be where the steady-state tem- 
perature is of the same order of the initial equilibrium 
temperature. In this case, y* and £* are not independent, 
and must satisfy a certain scaling relation. This relation 
is obtained by requiring y* 2 ~ 1 — e 2 , which leads to 

y*2 

c *^iV (20) 

We adopt this scaling relation, with T* = 10~ 4 , which is 
determined by the choice of k*. 

In Fig. 1, the results for the density time-correlator 
^> q (t) and the density-current time-correlator l P 9 (f) are 
shown. The shear rate is /To = 10~ 2 , the volume fraction 
<p is e = (<p — <p c )/<Pe = +10~ 3 , where (p c ~ 0.516 is the 



critical MCT transition point in equilibrium, and £ is var- 
ied, which corresponds to varying e. It can be seen that 
<& q (t) is almost coincident for < 1.0, which implies 
that the effect of dissipation is negligible. However, for 
> 1.0, which corresponds to the case of e < 0.9999, 
the plateau of ^> q (t) is destroyed, and & q (t) converges 
to the Debye-type relaxation curve at £* ~ 2.0, i.e. e ~ 
0.9998. This drastic demolishing of the plateau is mainly 
caused by the dissipative memory kernels Mq~\x) and 
M^ 3 '(t), Eqs. (17) and (18), which are of the opposite 

sign compared to Mq\x), as mentioned previously. In- 
tuitively, this feature can be regarded as the exhaustion 
of the "cage" due to inelastic collisions, which enables 
the particles to escape. Note that the demolishing of the 
plateau is accomplished at t ~ lOfo, which is far below 
the a-relaxation time x a ~ 10 2 To, where the effect of 
shear becomes significant. This may justify the applica- 
tion of the isotropic approximation. 

Furthermore, it is remarkable that, complementarily 
to ^> q (t), a plateau emerges in the tail of ^(f) for 
£* > 1.0. This implies that the density-current corre- 
lation is essential to the demolishing of the plateau of 
®q(t), which supports the role of the dissipative memory 

kernels, A/f^r) andAjfV)- 

The obtained result is compatible with the result 
from molecular dynamics, which also indicates that the 
plateau of the density time-correlator dissapears for e < 

0. 99, while it recovers for the nearly elastic case [29]. 

CONCLUDING REMARKS 

It is notable that the projection onto the density-current 
modes, Eq. (5), plays an essential role. This is in contrast 
to the sheared thermostatted systems, where the effect of 
^ n j(t) has been proved to be negligible [26]. This dif- 
ference resides in the nature of the dissipative interaction, 

1. e. whether it is single-body or two-body. 

According to the integration-through-transient scheme 
of Ref. [30], it is able to derive a formula for the steady- 
state shear stress in terms of time-correlators. Then, we 
can obtain the constitutive equation for the shear stress, 
which is to be compared to the result of molecular dy- 
namics [6, 7]. This issue, with an analysis near the jam- 
ming transition point, will be published elsewhere [25]. 

In conclusion, we have successfully formulated a 
MCT for sheared granular liquids. We have demonstrated 
that the plateau of the density time-correlator disappears 
for e < 0.9999, as observed in molecular dynamics sim- 
ulations. This destruction of the plateau results from the 
effect of the density-current correlation. 
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